Per Baubec request (October 17th 2019). Largely based upon cph_quickchanges_simpler.Rmd; older report available at taupo

Requires extract_motifs_frequency_from_bam_binary.sh to be run first, as wrapped up by (no coverage filter) motif_extract_run_nov_2019_no_coverage_filtering.sh.

Motifs are classified as methylated or unmethylated and counted, without aggregating average methylation values per sample to avoid batch/genotype effects.

Config

Data genotypes/annotations

Please check/update annotations:

DT::datatable(as.data.frame(annot), 
              extensions = c("Buttons", "FixedColumns"),
              rownames = FALSE, 
              options = list(dom = "Bfrtip",
                             scrollX = TRUE, 
                             fixedColumns = list(leftColumns = 1),
                             buttons = c("csv", "excel")))

Note that the merged datasets artificially gather samples from different runs, as follows:

DT::datatable(as.data.frame(tmp), 
              extensions = c("Buttons", "FixedColumns"),
              rownames = FALSE, 
              options = list(dom = "Bfrtip",
                             scrollX = TRUE, 
                             fixedColumns = list(leftColumns = 1),
                             buttons = c("csv", "excel")))

Overview

Aim

  • Hypothesis: DNMTs (DNA meth transferases) show sequence specificity
  • Test: check whether methylated/unmethylated motifs representation

Design

knitr::include_graphics("/home/imallona/src/cg_context/stranded_dnameth.png")

  • Flow
    • Call DNA methylation for every CpG in a strand-specific manner
    • Collapse counts (methylated/unmethylated) for each context 6-mer NNNCGNNN
    • Get nucleotide proportions for each position (i.e. 1 to 8) for methylated and unmethylated motifs separately
    • Get a log2 ratio

Datasets load

sources <- grep('cg_context|neuro', list.files(WD, recursive = TRUE, pattern = 'motif_counts_*'),
                invert = TRUE, value = TRUE)

Removing the merged 3a1, that only contains a single run

sources <- grep('d3a1_merged', sources, value = TRUE, invert = TRUE)
## dictionary of the whole set of motifs
nucl <- c('T', 'C', 'G', 'A')

mdict <- list('cg' = NULL, 'ch' = NULL)
mdict$cg <- apply(expand.grid('1' = nucl, '2' = nucl, '3' = nucl,
                               '4' = 'C',  '5' = 'G',
                               '6' = nucl, '7' = nucl, '8' = nucl ),
                   1,
                   function(x) paste(x, collapse = ''))



mdict$ch <- apply(expand.grid('1' = nucl, '2' = nucl, '3' = nucl,
                               '4' = 'C',  '5' = c('A', 'T', 'C'),
                               '6' = nucl, '7' = nucl, '8' = nucl),
                   1,
                   function(x) paste(x, collapse = ''))

for (item in names(mdict)) {
    fd <- mdict[[item]]
    fd <- sort(as.character(fd))
    fd <- data.frame(motif = fd, count = 0)
    rownames(fd) <- fd$motif
    mdict[[item]] <- fd
    rm(fd)
}

## filling with zeroes the motifs that are not present in some datasets

for (context in c('cg', 'ch')) {
    for (genotype in dirname(sources)) {
        for (state in c('meth', 'unmeth')) {
            curr <- d[[context]][[genotype]][[state]]

            mis <- setdiff(rownames(mdict[[context]]), rownames(curr))
            curr <- rbind(curr, mdict[[context]][mis,])
            ## lexicographical sorting
            curr <- curr[order(as.character(curr$motif)),]
            d[[context]][[genotype]][[state]] <- curr
            
            
            stopifnot(nrow(curr) == nrow(mdict[[context]]))
            rm(curr, mis)
        }
    }
}

Saving the data as an R object with the following structure:

  • context (cg or ch)
  • sample (i.e. qko)
  • status (meth or unmeth)
  • dataframe the columns motif, number of instances

Download

save(x = d, file = file.path(WD, 'counts_nested_list.RData'))
## str(d)

Seqlogo-based analysis

Based in nucleotide deviations from the 0.25:0.25:0.25:0.25 expectations for each nucleotide for each position.

proportion <- function(x){
   rs <- sum(x);
   return(x / rs);
}


deviation_from_proportion <- function(x){
   rs <- sum(x);
   ## return((x / rs) - 0.25)
   return((x / rs) / 0.25)
}


log2ratio_proportion <- function(x){
   rs <- sum(x);
   ## return((x / rs) - 0.25)
   return(log2((x / rs) - 0.25))
}

Data processing

Nucleotide weights are computed as if a PWM. We first get the nucleotide representation in methylated xor unmethylated motifs, and them compute the log2 enrichment of the meth/unmeth ratio to get a positive weight for overrepresentations. A step by step analysis is provided below (header Processing steps in detail).

outer <- list()
for (context in names(d)) {
    melted <- list()
    for (genotype in names(d[[context]])) {
        ## print(genotype)
        
        selected <- d[[context]][[genotype]]$meth

        for_pwm <- data.frame(A = rep(NA, nchar(selected$motif[1])),
                              C = NA,
                              G = NA,
                              T = NA)

        for (i in 1:nchar(selected$motif[1])) {
            id <- paste0('d',i)
            selected[, id] <- substr(selected$motif, i, i)
            for (nt in c('A', 'T', 'C', 'G')) {
                tsum <- 0
                tsum <- sum(selected[selected[,id] == nt,'count'])
                for_pwm[i, nt] <- tsum 
                
            }
        }
        
        pwm <- makePWM(apply(for_pwm, 1, proportion))
        pwm_meth <- pwm

        ## same for unmeth
        selected <- d[[context]][[genotype]]$unmeth

        for_pwm <- data.frame(A = rep(NA, nchar(selected$motif[1])),
                              C = NA,
                              G = NA,
                              T = NA)

        for (i in 1:nchar(selected$motif[1])) {
            id <- paste0('d',i)
            selected[, id] <- substr(selected$motif, i, i)
            for (nt in c('A', 'T', 'C', 'G')) {
                tsum <- 0
                tsum <- sum(selected[selected[,id] == nt,'count'])
                for_pwm[i, nt] <- tsum                 
            }
        }
        
        pwm <- makePWM(apply(for_pwm, 1, proportion))   
        pwm_unmeth <- pwm

        deviations <- log2(pwm_meth@pwm /pwm_unmeth@pwm)
        data.m <- melt(deviations)
        data.m$genotype <- genotype
        melted[[genotype]] <- data.m
        outer[[context]] <- melted        
    }
}

Processing steps in detail

Let’s visualize the calculations for sample tko+d3a2_merged and a context CG

context <- 'cg' ; genotype <- 'tko+d3a2_merged'

selected <- d[[context]][[genotype]]$meth

for_pwm <- data.frame(A = rep(NA, nchar(selected$motif[1])),
                      C = NA, G = NA,  T = NA)

for (i in 1:nchar(selected$motif[1])) {
    id <- paste0('d',i)
    selected[, id] <- substr(selected$motif, i, i)
    for (nt in c('A', 'T', 'C', 'G')) {
        tsum <- 0
        tsum <- sum(selected[selected[,id] == nt,'count'])
        for_pwm[i, nt] <- tsum 
        
    }
}

pwm <- makePWM(apply(for_pwm, 1, proportion))

cat(sprintf('First we calculate the nucleotide proportions for each nucleotide in methylated motifs\n'))
## First we calculate the nucleotide proportions for each nucleotide in methylated motifs
print(kable(pwm@pwm))
## 
## 
##               1           2           3    4    5           6           7           8
## ---  ----------  ----------  ----------  ---  ---  ----------  ----------  ----------
## A     0.2653976   0.1846908   0.3127135    0    0   0.2078319   0.2464551   0.3020369
## C     0.2570343   0.3280534   0.2994549    1    0   0.2979749   0.2892987   0.2562254
## G     0.2350310   0.1723616   0.1820777    0    1   0.1975226   0.2248220   0.1793862
## T     0.2425372   0.3148942   0.2057539    0    0   0.2966707   0.2394243   0.2623514
cat("\n")
pwm_meth <- pwm

## same for unmeth
selected <- d[[context]][[genotype]]$unmeth

for_pwm <- data.frame(A = rep(NA, nchar(selected$motif[1])),
                      C = NA,
                      G = NA,
                      T = NA)

for (i in 1:nchar(selected$motif[1])) {
    id <- paste0('d',i)
    selected[, id] <- substr(selected$motif, i, i)
    for (nt in c('A', 'T', 'C', 'G')) {
        tsum <- 0
        tsum <- sum(selected[selected[,id] == nt,'count'])
        for_pwm[i, nt] <- tsum                 
    }
}
pwm <- makePWM(apply(for_pwm, 1, proportion))
        
cat('We do the same for unmethylated motifs\n')
## We do the same for unmethylated motifs
print(kable(pwm@pwm))
## 
## 
##               1           2           3    4    5           6           7           8
## ---  ----------  ----------  ----------  ---  ---  ----------  ----------  ----------
## A     0.2777496   0.2696658   0.2702918    0    0   0.2470913   0.2426310   0.2340711
## C     0.2483383   0.2674804   0.2567395    1    0   0.1964955   0.2063492   0.2191563
## G     0.2354485   0.2402713   0.2302734    0    1   0.2880217   0.3020977   0.2760242
## T     0.2384636   0.2225825   0.2426953    0    0   0.2683915   0.2489221   0.2707485
cat("\n")      
pwm_unmeth <- pwm

cat(sprintf('We finally calculate the log2 (meth/unmeth) ratios\n'))
## We finally calculate the log2 (meth/unmeth) ratios
deviations <- log2(pwm_meth@pwm /pwm_unmeth@pwm)

print(kable(deviations))    
## 
## 
##                1            2            3     4     5            6            7            8
## ---  -----------  -----------  -----------  ----  ----  -----------  -----------  -----------
## A     -0.0656295   -0.5460601    0.2103238   NaN   NaN   -0.2496268    0.0225606    0.3677781
## C      0.0496542    0.2944975    0.2220335     0   NaN    0.6006940    0.4874724    0.2254537
## G     -0.0025610   -0.4792261   -0.3387937   NaN     0   -0.5441600   -0.4262323   -0.6217255
## T      0.0244368    0.5005267   -0.2382263   NaN   NaN    0.1445234   -0.0561250   -0.0454525

Coloring by annotations

All experiments

for (context in names(outer)) {
    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')
    melted2 <- merge(melted2, annot, by = 'sample', all.x = TRUE)
   
    print(ggplot(melted2, aes(position, deviation, fill=as.factor(sample)))+
          geom_bar(position="dodge",stat="identity") +
          scale_x_continuous(breaks = 1:8) +
          facet_wrap(~id+nucleotide, ncol= 4) +
          xlab("position") +          
          ylab('log2fc meth/unmeth') +
          theme(legend.position = 'none') + ggtitle(context))
}
## Warning: Removed 228 rows containing missing values (geom_bar).

## Warning: Removed 152 rows containing missing values (geom_bar).

Per sample quality checks

Coverage (total number of methylated plus unmethylated motifs detected).

counts <- list()
for (context in names(d)) {
    tmp <- list()
    for (sample in names(d[[context]])) {
        meth <- d[[context]][[sample]]$meth
        meth$sample <- sample
        meth$context <- context
        meth$status <- 'meth'
        
        unmeth <- d[[context]][[sample]]$unmeth

        unmeth$sample <- sample
        unmeth$context <- context
        unmeth$status <- 'unmeth'
        
        tmp[[sample]] <- rbind(meth, unmeth)
       
        
    }

    counts[[context]] <- do.call(rbind.data.frame, tmp)
  
}
for (context in names(counts)) {
    sums <- aggregate(counts[[context]]$count,
                      by=list(counts[[context]]$sample,
                              counts[[context]]$status), sum)
    colnames(sums) <- c('sample', 'status', 'count')

    sums <- merge(sums, annot, by ='sample')
    print(ggplot(sums,
           aes(x = sample, y = count , fill = factor(id)))+
        geom_bar(stat="identity",position="dodge") +
        scale_y_continuous(trans='log10') +
        facet_wrap(~status , nrow = 1) + 
        xlab("Sample") + ylab("Cum. number of motifs detected") +
        theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
        ggtitle(context))
}

And by sequencing technology

for (context in names(counts)) {
    sums <- aggregate(counts[[context]]$count,
                      by=list(counts[[context]]$sample,
                              counts[[context]]$status), sum)
    colnames(sums) <- c('sample', 'status', 'count')

    sums <- merge(sums, annot, by ='sample')
    print(ggplot(sums,
           aes(x = sample, y = count , fill = factor(id)))+
        geom_bar(stat="identity",position="dodge") +
        scale_y_continuous(trans='log10') +
        facet_wrap(~status , nrow = 1) + 
        xlab("Sample") + ylab("Cum. number of motifs detected") +
        theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
        ggtitle(context))
}

Total number of motifs detected at least once (methylated or unmethylated).

Note that for CG we expect 4^6=4096, and for CH 4^6*3=12288.

for (context in names(counts)) {
    detects <- aggregate(counts[[context]]$count,
                      by=list(counts[[context]]$sample,
                              counts[[context]]$status), function(x) sum(x >=1))
    colnames(detects) <- c('sample', 'status', 'count')

    detects <- merge(detects, annot, by ='sample')
    print(ggplot(detects,
           aes(x = sample, y = count , fill = factor(id)))+
        geom_bar(stat="identity",position="dodge") +
        scale_y_continuous(trans='log10') +
        facet_wrap(~status , nrow = 1) + 
        xlab("Sample") + ylab("Number of motifs detected (at least once)") +
        theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
        ggtitle(context))
}

Per nucleotide and sample

for (context in names(d)) {
    
    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')
    
    print(ggplot(melted2, aes(position, deviation, fill=as.factor(sample)))+
          geom_bar(position="dodge",stat="identity") +
          scale_x_continuous(breaks = 1:8) +
          facet_wrap(~nucleotide,nrow = 2) + xlab("position") +
          ylab('log2fc meth/unmeth') +
          ggtitle(context) +
          guides(fill = guide_legend(nrow = length(unique(melted2$sample)))) + 
          theme(legend.position="bottom", legend.title=element_text(size=5), 
                legend.text=element_text(size=8), legend.key.size = unit(0.7, "line")))
    
}
## Warning: Removed 228 rows containing missing values (geom_bar).

## Warning: Removed 152 rows containing missing values (geom_bar).

for (context in names(d)) {
    idx <- names(sort(sapply(outer[[context]], function(x)
        return(x[x$Var1 == 'C' & x$Var2 == 2, 'value']))))

    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')
    melted2$sample <- factor(x = melted2$sample, levels = idx)

    print(ggplot(melted2, aes(position, deviation, fill=as.factor(sample)))+
          geom_bar(position="dodge",stat="identity") +
          scale_x_continuous(breaks = 1:8) +
          facet_wrap(~nucleotide,nrow = 2) + xlab("position") +
          ylab('log2fc meth/unmeth') +
          ggtitle(context) +
          guides(fill = guide_legend(nrow = length(unique(melted2$sample)))) + 
          theme(legend.position="bottom", legend.title=element_text(size=5), 
                legend.text=element_text(size=8), legend.key.size = unit(0.7, "line")))
    
   
}

Per nucleotide, celltype and annotation

for (context in names(d)) {
    tmp <- as.numeric(annot[names(outer$cg),'id_short'])
    names(tmp) <- names(outer$cg)

    ## actually this is not being used
    idx <- names(sort(tmp))

    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')

    annot$sample <- rownames(annot)
    melted2 <- merge(melted2, annot, by = 'sample')
    
    melted2$sample <- factor(x = melted2$sample, levels = idx)

    print(ggplot(melted2, aes(position, deviation, fill=as.factor(id_short)))+
          geom_bar(position="dodge",stat="identity") +
          scale_x_continuous(breaks = 1:8) +
          facet_wrap(~celltype + nucleotide, nrow = 2) + xlab("position") +
          ylab('log2fc meth/unmeth') +
          ggtitle(context))
    
    
}
## Warning: Removed 228 rows containing missing values (geom_bar).

## Warning: Removed 152 rows containing missing values (geom_bar).

Per celltype, annotation and nucleotide

for (context in names(d)) {
   
    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')

    annot$sample <- rownames(annot)
    melted2 <- merge(melted2, annot, by = 'sample')
    
    melted2$sample <- factor(x = melted2$sample, levels = idx)

    print(ggplot(melted2, aes(nucleotide, deviation, fill=as.factor(position)))+
          geom_bar(position="dodge",stat="identity") +
          facet_wrap(~celltype + id_short, nrow = 2) + xlab("genotype") +
          ylab('log2fc meth/unmeth') +
          ggtitle(context))


}
## Warning: Removed 228 rows containing missing values (geom_bar).

## Warning: Removed 152 rows containing missing values (geom_bar).

for (context in names(d)) {
   
    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')

    annot$sample <- rownames(annot)
    melted2 <- merge(melted2, annot, by = 'sample')
    
    melted2$sample <- factor(x = melted2$sample, levels = idx)

    print(ggplot(melted2, aes(nucleotide, deviation, fill=as.factor(position)))+
          geom_bar(position="dodge",stat="identity") +
          facet_wrap(~celltype + id_short, nrow = 2) + xlab("id") +
          ylab('log2fc meth/unmeth') +
          ggtitle(context))

}
## Warning: Removed 228 rows containing missing values (geom_bar).

## Warning: Removed 152 rows containing missing values (geom_bar).

Per sample and nucleotide

for (context in names(d)) {
    
    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')

    annot$sample <- rownames(annot)
    melted2 <- merge(melted2, annot, by = 'sample')
    
    melted2$sample <- factor(x = melted2$sample, levels = idx)

        
    print(ggplot(melted2, aes(position, deviation, fill=as.factor(nucleotide)))+
          geom_bar(position="dodge",stat="identity") +
          scale_x_continuous(breaks = 1:8) +
          facet_wrap(~sample, nrow=12) + xlab("position") + ylab('log2fc meth/unmeth') +
          ggtitle(context))

    
}
## Warning: Removed 228 rows containing missing values (geom_bar).

## Warning: Removed 152 rows containing missing values (geom_bar).

Naive Bayes

(ongoing)

Subset

Plotting some genotypes only

## Warning: Removed 66 rows containing missing values (geom_bar).

## Warning: Removed 66 rows containing missing values (geom_bar).

## Warning: Removed 44 rows containing missing values (geom_bar).

## Warning: Removed 44 rows containing missing values (geom_bar).

Clustering

Simple approach: cluster the overall rate of meth/unmeth per motif (this is biased by methylation level).

ratios <- sapply(d$cg, function(x) return(x$meth$count))/
    (sapply(d$cg, function(x) return(x$unmeth$count)) + sapply(d$cg, function(x) return(x$meth$count)) ) 

rownames(ratios) <- d$cg[[1]]$unmeth$motif

ratios <- ratios[,c(targets, grep('*merged*', colnames(ratios), value = TRUE))]
    
hc <- hclust(dist(t(ratios)), method = 'ward.D')
plot(hc)

Annotated clustering

## hc2 <- hc

## hc2$labels <- sprintf('%s\n(%s)', hc$labels, annot[hc$labels,'id'])
dend <- as.dendrogram(hc)

## labels_colors(dend) <- as.numeric(as.factor( annot[hc$labels,'id']))
par(mar = c(4,10,1,22), xpd=TRUE)
plot(dend, horiz = TRUE)
colored_bars(cbind(as.numeric(as.factor(annot[hc$labels,c('id_short')])),
                   as.numeric(as.factor(annot[hc$labels,c('tech')])),
                   as.numeric(as.factor(annot[hc$labels,c('celltype')]))), dend, horiz = TRUE,
             rowLabels = c('id', 'tech','celltype'))

legend("bottomleft",inset=c(-0.5,0), title ='id',
       legend = levels(as.factor(annot[hc$labels,'id_short'])),
       fill = 1:nlevels(as.factor(annot[hc$labels,'id'])))

legend("topleft", inset=c(-0.5,0), title = 'tech',
       legend = levels(as.factor(annot[hc$labels,'tech'])),
       fill = 1:nlevels(as.factor(annot[hc$labels,'tech'])))

rownames(annot) <- annot$sample
annot <- annot[,-1]

pheatmap(ratios, cluster.rows = TRUE, cluster.cols = TRUE, annotation_col = annot,
         scale = 'none', show_rownames = FALSE)

## note the scaling and centering by column
pheatmap(ratios, cluster.rows = TRUE, cluster.cols = TRUE, annotation_col = annot,
         scale = 'column', show_rownames = FALSE)

Heatmap representations

Per Tuncay request 11 Nov 2019

for (context in names(d)) {
    
    melted2 <- do.call(rbind.data.frame, outer[[context]])
    colnames(melted2) <- c('nucleotide', 'position', 'deviation', 'sample')

    annot$sample <- rownames(annot)
    melted2 <- merge(melted2, annot, by = 'sample')
    
    melted2$sample <- factor(x = melted2$sample, levels = idx)

    for (sample in levels(melted2$sample)) {
        tmp <- melted2[melted2$sample == sample,]

        wide <- dcast(tmp[,c('deviation', 'nucleotide', 'position')],  nucleotide ~ position,
                      value.var = 'deviation')
        rownames(wide) <- wide$nucleotide
        wide <- wide[,-1]

        ## while (!is.null(dev.list())) dev.off()
        ## cat(sprintf('\n%s %s\n\n\n', context, sample))
        cat('\n\n\n\n')
        pheatmap(wide, cluster_rows = FALSE, cluster_cols = FALSE,
                       color = colorRampPalette(c('red', 'white', 'blue'))(60),
                       breaks = seq(-3, 3, 0.1), main = sprintf('%s\n%s', context, sample))
       
    }
}

Timestamp

date()
## [1] "Mon Nov 18 17:38:52 2019"
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/local/R/R-3.6.0/lib/libRblas.so
## LAPACK: /usr/local/R/R-3.6.0/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
##  [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
##  [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] dendextend_1.12.0 viridis_0.5.1     viridisLite_0.3.0
##  [4] naivebayes_0.9.6  DT_0.9            gridExtra_2.3    
##  [7] ggseqlogo_0.1     ggfortify_0.4.7   ggplot2_3.2.1    
## [10] seqLogo_1.50.0    pheatmap_1.0.12   reshape2_1.4.3   
## [13] knitr_1.25       
## 
## loaded via a namespace (and not attached):
##  [1] pkgload_1.0.2      tidyr_1.0.0        jsonlite_1.6      
##  [4] shiny_1.3.2        assertthat_0.2.1   highr_0.8         
##  [7] stats4_3.6.0       yaml_2.2.0         remotes_2.1.0     
## [10] sessioninfo_1.1.1  pillar_1.4.2       backports_1.1.4   
## [13] glue_1.3.1         digest_0.6.20      RColorBrewer_1.1-2
## [16] promises_1.0.1     colorspace_1.4-1   htmltools_0.3.6   
## [19] httpuv_1.5.2       plyr_1.8.4         pkgconfig_2.0.2   
## [22] devtools_2.2.0     purrr_0.3.2        xtable_1.8-4      
## [25] scales_1.0.0       processx_3.4.1     later_0.8.0       
## [28] tibble_2.1.3       usethis_1.5.1      ellipsis_0.2.0.1  
## [31] withr_2.1.2        lazyeval_0.2.2     cli_1.1.0         
## [34] magrittr_1.5       crayon_1.3.4       mime_0.7          
## [37] memoise_1.1.0      evaluate_0.14      ps_1.3.0          
## [40] fs_1.3.1           pkgbuild_1.0.5     tools_3.6.0       
## [43] prettyunits_1.0.2  lifecycle_0.1.0    stringr_1.4.0     
## [46] munsell_0.5.0      callr_3.3.1        compiler_3.6.0    
## [49] rlang_0.4.0        htmlwidgets_1.3    crosstalk_1.0.0   
## [52] labeling_0.3       rmarkdown_1.15     testthat_2.2.1    
## [55] gtable_0.3.0       codetools_0.2-16   R6_2.4.0          
## [58] dplyr_0.8.3        zeallot_0.1.0      rprojroot_1.3-2   
## [61] desc_1.2.0         stringi_1.4.3      Rcpp_1.0.2        
## [64] vctrs_0.2.0        png_0.1-7          tidyselect_0.2.5  
## [67] xfun_0.9
devtools::session_info()
## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.0 (2019-04-26)
##  os       Ubuntu 16.04.6 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language en_CA:en                    
##  collate  en_CA.UTF-8                 
##  ctype    en_CA.UTF-8                 
##  tz       Europe/Zurich               
##  date     2019-11-18                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package      * version date       lib source        
##  assertthat     0.2.1   2019-03-21 [1] CRAN (R 3.6.0)
##  backports      1.1.4   2019-04-10 [1] CRAN (R 3.6.0)
##  callr          3.3.1   2019-07-18 [1] CRAN (R 3.6.0)
##  cli            1.1.0   2019-03-19 [1] CRAN (R 3.6.0)
##  codetools      0.2-16  2018-12-24 [3] CRAN (R 3.6.0)
##  colorspace     1.4-1   2019-03-18 [1] CRAN (R 3.6.0)
##  crayon         1.3.4   2017-09-16 [1] CRAN (R 3.6.0)
##  crosstalk      1.0.0   2016-12-21 [1] CRAN (R 3.6.0)
##  dendextend   * 1.12.0  2019-05-11 [1] CRAN (R 3.6.0)
##  desc           1.2.0   2018-05-01 [1] CRAN (R 3.6.0)
##  devtools       2.2.0   2019-09-07 [1] CRAN (R 3.6.0)
##  digest         0.6.20  2019-07-04 [1] CRAN (R 3.6.0)
##  dplyr          0.8.3   2019-07-04 [1] CRAN (R 3.6.0)
##  DT           * 0.9     2019-09-17 [1] CRAN (R 3.6.0)
##  ellipsis       0.2.0.1 2019-07-02 [1] CRAN (R 3.6.0)
##  evaluate       0.14    2019-05-28 [1] CRAN (R 3.6.0)
##  fs             1.3.1   2019-05-06 [1] CRAN (R 3.6.0)
##  ggfortify    * 0.4.7   2019-05-26 [2] CRAN (R 3.6.0)
##  ggplot2      * 3.2.1   2019-08-10 [1] CRAN (R 3.6.0)
##  ggseqlogo    * 0.1     2017-07-25 [2] CRAN (R 3.6.0)
##  glue           1.3.1   2019-03-12 [1] CRAN (R 3.6.0)
##  gridExtra    * 2.3     2017-09-09 [1] CRAN (R 3.6.0)
##  gtable         0.3.0   2019-03-25 [1] CRAN (R 3.6.0)
##  highr          0.8     2019-03-20 [1] CRAN (R 3.6.0)
##  htmltools      0.3.6   2017-04-28 [1] CRAN (R 3.6.0)
##  htmlwidgets    1.3     2018-09-30 [1] CRAN (R 3.6.0)
##  httpuv         1.5.2   2019-09-11 [1] CRAN (R 3.6.0)
##  jsonlite       1.6     2018-12-07 [1] CRAN (R 3.6.0)
##  knitr        * 1.25    2019-09-18 [1] CRAN (R 3.6.0)
##  labeling       0.3     2014-08-23 [1] CRAN (R 3.6.0)
##  later          0.8.0   2019-02-11 [1] CRAN (R 3.6.0)
##  lazyeval       0.2.2   2019-03-15 [1] CRAN (R 3.6.0)
##  lifecycle      0.1.0   2019-08-01 [1] CRAN (R 3.6.0)
##  magrittr       1.5     2014-11-22 [1] CRAN (R 3.6.0)
##  memoise        1.1.0   2017-04-21 [1] CRAN (R 3.6.0)
##  mime           0.7     2019-06-11 [1] CRAN (R 3.6.0)
##  munsell        0.5.0   2018-06-12 [1] CRAN (R 3.6.0)
##  naivebayes   * 0.9.6   2019-06-03 [2] CRAN (R 3.6.0)
##  pheatmap     * 1.0.12  2019-01-04 [1] CRAN (R 3.6.0)
##  pillar         1.4.2   2019-06-29 [1] CRAN (R 3.6.0)
##  pkgbuild       1.0.5   2019-08-26 [1] CRAN (R 3.6.0)
##  pkgconfig      2.0.2   2018-08-16 [1] CRAN (R 3.6.0)
##  pkgload        1.0.2   2018-10-29 [1] CRAN (R 3.6.0)
##  plyr           1.8.4   2016-06-08 [1] CRAN (R 3.6.0)
##  png            0.1-7   2013-12-03 [1] CRAN (R 3.6.0)
##  prettyunits    1.0.2   2015-07-13 [1] CRAN (R 3.6.0)
##  processx       3.4.1   2019-07-18 [1] CRAN (R 3.6.0)
##  promises       1.0.1   2018-04-13 [1] CRAN (R 3.6.0)
##  ps             1.3.0   2018-12-21 [1] CRAN (R 3.6.0)
##  purrr          0.3.2   2019-03-15 [1] CRAN (R 3.6.0)
##  R6             2.4.0   2019-02-14 [1] CRAN (R 3.6.0)
##  RColorBrewer   1.1-2   2014-12-07 [1] CRAN (R 3.6.0)
##  Rcpp           1.0.2   2019-07-25 [1] CRAN (R 3.6.0)
##  remotes        2.1.0   2019-06-24 [1] CRAN (R 3.6.0)
##  reshape2     * 1.4.3   2017-12-11 [1] CRAN (R 3.6.0)
##  rlang          0.4.0   2019-06-25 [1] CRAN (R 3.6.0)
##  rmarkdown      1.15    2019-08-21 [1] CRAN (R 3.6.0)
##  rprojroot      1.3-2   2018-01-03 [1] CRAN (R 3.6.0)
##  scales         1.0.0   2018-08-09 [1] CRAN (R 3.6.0)
##  seqLogo      * 1.50.0  2019-05-02 [2] Bioconductor  
##  sessioninfo    1.1.1   2018-11-05 [1] CRAN (R 3.6.0)
##  shiny          1.3.2   2019-04-22 [1] CRAN (R 3.6.0)
##  stringi        1.4.3   2019-03-12 [1] CRAN (R 3.6.0)
##  stringr        1.4.0   2019-02-10 [1] CRAN (R 3.6.0)
##  testthat       2.2.1   2019-07-25 [1] CRAN (R 3.6.0)
##  tibble         2.1.3   2019-06-06 [1] CRAN (R 3.6.0)
##  tidyr          1.0.0   2019-09-11 [1] CRAN (R 3.6.0)
##  tidyselect     0.2.5   2018-10-11 [1] CRAN (R 3.6.0)
##  usethis        1.5.1   2019-07-04 [1] CRAN (R 3.6.0)
##  vctrs          0.2.0   2019-07-05 [1] CRAN (R 3.6.0)
##  viridis      * 0.5.1   2018-03-29 [1] CRAN (R 3.6.0)
##  viridisLite  * 0.3.0   2018-02-01 [1] CRAN (R 3.6.0)
##  withr          2.1.2   2018-03-15 [1] CRAN (R 3.6.0)
##  xfun           0.9     2019-08-21 [1] CRAN (R 3.6.0)
##  xtable         1.8-4   2019-04-21 [1] CRAN (R 3.6.0)
##  yaml           2.2.0   2018-07-25 [1] CRAN (R 3.6.0)
##  zeallot        0.1.0   2018-01-28 [1] CRAN (R 3.6.0)
## 
## [1] /home/Shared/Rlib/release-3.9-lib
## [2] /home/imallona/R/x86_64-pc-linux-gnu-library/3.6
## [3] /usr/local/R/R-3.6.0/library